# THIS SCRIPT REPLICATES "Data and Methods" AND "Results" SECTIONS
data <- read.csv("data.csv") 

cat("Table 2: Summary statistics", '\n')
t <- c(3, 4, 5, 11, 13, 15, 16)
cols <- paste("nelda", t, sep="")
names <- c("Opposition not allowed", "Only one party legal", "No choice on the ballot", "Concerns elections will not be free and fair", "Opponents prevented from running", "Opposition harassed", "Allegations of media bias")
tab <- cbind(paste(1:7, ".", sep = ""), paste("n", t, sep = ""), names = names ,roundr(apply(data[,  cols], 2, mean, na.rm=TRUE), 2), roundr(apply(data[,  cols], 2, sd, na.rm=TRUE), 2))
colnames(tab) <- c("Indicator", "Label", "Description", "Mean", "St. Dev.")
             
print(xtable(tab,
      caption = "Indicators of electoral manipulation.",
      label   = "tab:manipulation",
      align   = c("l", "c", "l", "l", "c", "c")),
      sanitize.text.function = identity,
      booktabs = TRUE,
      table.placement = 't!',
      include.rownames = FALSE,
      include.colnames = TRUE,
      file = "../tables/table2.tex"
)

cat("Figure 3: Identifying crisis threshold: will take a while", '\n')

Model <- function(x) paste("nelda24 ~ (1|country) + lambda + I(s18f1_lag > 0) + mirt + (growthc_t1 < ", as.character(x), ")", sep="")

fit <- function(x) glmer(Model(x), data = data, family = binomial, na.action = na.exclude)
x   <- seq(-2, 0, length = 501)
s   <- lapply(x, fit)

pdf("../graphs/gthreshold.pdf", width=10, height=6)
par(mar = c(5, 5, 1/2, 1/2), mgp = c(3.5,.75,0))
S <- unlist(lapply(s, logLik))
plot(x, S, type="l", lwd=2, xlab="Economic crisis threshold", ylab="Log-likelihood", cex.lab=1.5, las = 1)
fit2 <- glmer(nelda24 ~ (1|country) + lambda + I(s18f1_lag > 0) + mirt + growth_t1, data=data, family=binomial)
t <- grep("growthc_t1", names(fixef(s[[1]])))
p <- x[which(sapply(s, function(s) coefficients(summary(s))[t, 4] < 0.05))]
segments(p, min(S)-.05, p, max(S)+.05, col = "grey90", lwd=2)
abline(h=logLik(fit2), lty=3, lwd=3)
lines(x, unlist(lapply(s, logLik)), type="l", lwd=2.5)
text(-0.3, logLik(fit2) + .75, labels = "Log-likelihood of the model \n with a linear growth term", cex=1.25)
arrows(-.1, logLik(fit2) + .5, -.1, logLik(fit2), col="black", lwd=2)
dev.off()

cat("Figure 4: Reelection probabilities at s = 1", '\n')
fit <- glmer(Model(-1), data=data, family = binomial)
t <- c("growthc_t1 < -1TRUE", "lambda", "I(s18f1_lag > 0)TRUE", "mirt")
ci <- exp(confint(fit)[t, ])
m <- exp(fixef(fit))[t]

pdata <- data.frame(var = c("Crisis", "Coup threat", "Protests", "Manipulation"), Coefficient = m, cl = ci[, 1], cu = ci[,2])
pdata$var <- factor(pdata$var, c("Crisis", "Coup threat", "Protests", "Manipulation"))

ggplot(data = pdata, aes(x = var, y = Coefficient, ymin = cl, ymax = cu)) + geom_point(position = position_dodge(width = 0.2), size=4) + geom_errorbar(position = position_dodge(width = 0.2), width = 0.2, size=1) + scale_colour_manual(values = c("blue", "red")) + theme_classic() + geom_hline(yintercept=c(1), linetype="dotted") +  theme(axis.title.y=element_blank(), axis.text.x=element_text(size=15), axis.title.x=element_text(size=16, vjust = -.5), axis.text.y=element_text(size=20)) + ylab("Predicted relative odds with 95% CI's") + coord_flip() + geom_text(aes(y = Coefficient, x = (1:4)+.25, label = roundr(Coefficient, 2)))
ggsave("../graphs/losingodds.pdf", width=10, height=4)

cat("Table 3: Main results", '\n')
Model <- mirt ~ I(growthc_t1 < -1) + lambda + I(s18f1_lag > 0) + gdp.gle_t1 + agedem_t1 + nelda19 + legelec + factor(year) + (1|country)

fit1 <- lmer(Model, data = data, REML = FALSE)
fit2 <- update(fit1, . ~ . + partyt1 + militaryt1 + personalt1 + monarchyt1)
fit3 <- update(fit2, . ~ . + terror + (1|leader))

# Compute Hausman p-value: add manually to the above table
hausman.p <- function(fit) phtest_glmer(Mod1 = fit, Mod2 = plm(update(formula(fit), . ~ . - (1|country)), model = "within", index = "country", data = data))$p.value
p <- unlist(lapply(list(fit1, fit2), hausman.p))
hausman.p <- function(fit) phtest_glmer(Mod1 = fit, Mod2 = update(fit, . ~ . - (1|country) + factor(country)))$p.value
p <- roundr(c(p, hausman.p(fit3)), 2)


var.labels <- c("Intercept", "Crisis", "Coup threat", "Protests", "GDP per capita", "Democracy age", "Foreign aid", "Legislative elections", "Party institutions", "Military institutions", "Personalist institutions", "Monarchic institutions", "Political terror")


intable <- function(p, name) { 
      k <- paste("\\multicolumn{1}{r}{", p, "}", sep = "")
      paste(name, "&", paste(k, collapse = "&&"), "\\\\\n", collapse="")
}


add.line <- list(command = c(intable(p, "Hausman p-value"), intable(c("118", "112", "109/232"), "Countries/Leaders"), intable(rep("Yes", 3), "Country-effects"),  intable(rep("Yes", 3), "Year-effects"),intable(c("", "", "Yes"), "Leader-Effects")))

outtex(list(fit1, fit2, fit3), omit = c("year"), label = "tab:reg1", caption = "Regression estimates for electoral manipulation scores.", var.labels = var.labels, ci = list(method = "profile"), add.line = list(pos = rep(14, 5), command = add.line$command), below = FALSE, size = "small", space = "50pt", file = "../tables/table3.tex")

# Bayes factor comparison with polynomial-year model
bfact <- function(fit) BIC(fit) - BIC(update(fit, . ~ . - factor(year) + poly(year, 3)))
unlist(lapply(list(fit1, fit2, fit3), bfact))

cat("Figure 6: interaction with availability of polls", '\n')
fit.inter <- lmer(mirt ~ (1|country) + nelda19 + legelec + gdp.gle_t1 + agedem_t1 + ((growthc_t1 < -1) + lambda + (s18f1_lag > 0)):(nelda25==1) + (nelda25==1) + factor(year), data = data, REML = FALSE)

b   <- fixef(fit.inter)[which(grepl("growthc|lambda|s18f1", names(fixef(fit.inter))))]
cis <- confint(fit.inter)
cis <- cis[which(grepl("growthc|lambda|s18f1", rownames(cis))),]

pdata <- data.frame(var = rep(c("Crisis", "Coup threat", "Protest"), each=2), Polls=rep(c("Unavailable", "Available"), 3), Coefficient = b, cl = cis[, 1], cu=cis[, 2])

ggplot(data = pdata, aes(x = var, y = Coefficient, ymin = cl, ymax = cu, colour = Polls, shape = Polls)) + geom_point(position = position_dodge(width = 0.2), size=4) + geom_errorbar(position = position_dodge(width = 0.2), width = 0.2, size=1) + coord_flip() +  scale_colour_manual(values = c("grey50", "black")) + theme_classic() + geom_hline(xintercept=c(-3,0), linetype="dotted") +  theme(axis.title.y=element_blank(), axis.text.x=element_text(size=12), axis.title.x=element_text(size=14), axis.text.y=element_text(size=14), legend.key.size = unit(3, "cm"), legend.title=element_text(size=14), legend.text=element_text(size=14)) + scale_y_continuous(breaks=seq(-.4,.2,by = 0.1)) + ylab("Coefficient and 95% CI's")
ggsave("../graphs/polls.pdf", width=8, height=4)

cat("Appendix Table 1: with country fixed effects", "\n")
refit1 <- lm(update(formula(fit1), ~ . - (1|country) + country), data = data)
refit2 <- lm(update(formula(fit2), ~ . - (1|country) + country), data = data)
refit3 <- update(fit3, . ~ .  - (1|country) + country)

outtex(list(refit1, refit2, refit3), omit = c("year", "country"), caption = "Regression estimates for electoral manipulation scores with country and year fixed effects.", var.labels = var.labels, ci = list(method = "boot"), below = FALSE, label = "appendix_table1", table.placement = 'h!', size = "small", file = "../tables/appendix_table1.tex")
